/* perform complete mass function analysis
------------------------------------------
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "mfnanalyse.h"

void print_data(double q, double ftNL) {
 double M = 1.61e14;
 double m3r1=M3R(M, calA(0.),1.,1.);
 double m3r=M3R(M, calA(q*ftNL), q, ftNL);
 double rat=M3Rf(M, calA(q*ftNL), q, ftNL)/M3Rh(M, calA(q*ftNL), q, ftNL);
 printf("%g %g %g\n", m3r/m3r1, m3r, rat);
}

int main(int argc, char *argv[]) {
  //printf("USAGE: mfn q ftNL name FEEDER\n");
  double m3r1;
 m3r1 = M3R(1.61e14, calA(0), 1., 1.);
  printf("a "); print_data(0.1, 50000.);
  printf("cc "); print_data(0.00005, 100000000.);
  printf("fnl500 "); print_data(1., 500.);
  printf("aa "); print_data(0.11925, 20620.);
  printf("ccc "); print_data(0.00003, 100000000.);
  printf("eee "); print_data(0.00003, 80000000.);
  printf("fnl99 "); print_data(1., 99.);
  printf("sss "); print_data(0.00003, 65000000.);
  //return 0;

  printf("Some checks:\n-----------\n");
  printf("%g, %g, %g\n", d2R(1.61e14, 1, 1, 1), nuc(1.6, 1.61e14, 1, 1, 99, 0), M3Rf(1.61e14, calA(5000), 1, 5000));
  printf("%g, %g\n", ratio3f(1.61e14, calA(2000), 1.6, 1, 2000, 0), ratio4f(1.61e14, calA(0), 1.6, 1, 0, 0));
 //  return 0;
/*
  double aM=1.61e14;
  printf("3h: %g, 3f: %g, 4h: %g, 4f: %g, 5h: %g, 5f: %f\n", ratio3h(aM, AC, 1.6, 1, 99, 1), ratio3f(aM, AC, 1.6, 1, 500, 0), ratio4h(aM, AC, 1.6, 1, 99, 1), ratio4f(aM, AC, 1.6, 1, 500, 0), ratio5h(aM, AC, 1.6, 1, 99, 1), ratio5f(aM, AC, 1.6, 1, 500, 0));
  printf("I41: %g, I41D: %g\n", I4R1(aM), I4R1D(aM));
*/
  double chisq, chisqmin=5000., dcmin=10.;
  double dc, dcerror, numax;

  /* Results from PART 1 */
  /* PART 3: generate theory and error curves */
  FILE *fp;
  char fname[150], Opbase[200];
 // sprintf(Opbase, "/home1/02539/sza5154/feedersim/Theory/Output1024");
  sprintf(Opbase, "/home/sza5154/Research/feedersim/Theory/Output1024"); // this may need modification
  double M;
  int z;
  numax=12.;
 // find the best dc
/***
 for (dc=1.4; dc<1.5; dc=dc+0.01) {
    chisq = fit_mf(dc, 1, 99, "fnl99", numax)+fit_mf(dc, 1., 500., "fnl500", numax)+fit_mf(dc, 0.00003, 65000000., "sss", numax)+fit_mf(dc, 0.00003, 80000000., "eee", numax) + fit_mf(dc, 0.00003, 100000000, "ccc", numax);
    if (chisq<chisqmin) {
      chisqmin=chisq;
      dcmin=dc;
    }
      printf("dc: %g, chisq: %g\n", dcmin, chisqmin);    
  }
  
  printf("dcmin: %g\n", dcmin);
  return 0;
  
***/  
  
  dc=1.46;
  
  numax=12.;
  double totalchisq=0.;
  totalchisq+=fit_mf(dc, 1., 99., "fnl99", numax);
  totalchisq+=fit_mf(dc, 1., 500., "fnl500", numax);
  totalchisq+=fit_mf(dc, 0.00003, 65000000., "sss",numax);
  totalchisq+=fit_mf(dc, 0.00003, 80000000., "eee", numax);
  totalchisq+=fit_mf(dc, 0.00003, 100000000, "ccc", numax);
  totalchisq+=fit_mf(dc, 0.00005, 100000000, "cc", numax);
//  totalchisq+=fit_mf(dc, 0.1, 50000,  "a", numax);
//  totalchisq+=fit_mf(dc, 0.119246, 20620, "aa", numax);
  
  printf("total chisq: %g\n", totalchisq);
  generate_theory_curves(dc, 0.1, 50000, Opbase, "a", 0);
  generate_theory_curves(dc, 0.119246, 20620, Opbase, "aa", 0);
  
//  generate_theory_m3(dc,  Opbase);
/***
  dc=1.686;
  generate_theory_curves(dc, 1., 99., Opbase, "fnl99",  1);
  generate_theory_curves(dc, 1., 500., Opbase,  "fnl500", 1);
  generate_theory_curves(dc, 0.00003, 65000000.,  Opbase, "sss", 1);
  generate_theory_curves(dc, 0.00003, 80000000., Opbase, "eee", 1);
  generate_theory_curves(dc, 0.00003, 100000000., Opbase, "ccc", 1);
  generate_theory_curves(dc, 0.00005, 100000000., Opbase, "cc", 1);
***/
}

double fit_mf(double dc, double q, double ftNL, char* name, double numax) {
  double chisq, chisqmin=5000., dcmin=10.;
 double  dcerror;

  /* PART 0: find the best fit dc using nu<3.5 from H99, H500, eee, and ccc. */

  numax=12.;
  double cf, cfmin, cferror;
  double cmin, cmax;
  if (name == "cc" || name =="a") {
    cmin=1.1; cmax=1.3;
  } else if (name=="fnl99" || name=="eee") {
    cmin=1.0; cmax=1.05;
  } else {
    cmin=1.0; cmax=1.1;
  }
for (cf=cmin; cf<cmax; cf=cf+0.001) {
//  for (dc=1.20;dc<1.80; dc=dc+0.01) {
      chisq = mfn(dc, q, ftNL, name, numax, cf);
//      printf("%s, dc: %g, cf: %g chisq: %g\n", name, dc,cf, chisq);
      if (chisq<chisqmin) {
        chisqmin=chisq;
        cfmin=cf;
   }
}

  cf=cfmin;
  chisq=chisqmin;
    while (chisq<chisqmin+1.0) {
    chisq=mfn(dc, q, ftNL, name, numax, cf);
    cferror=cf-cfmin;
    cf=cf+0.001;
  }
  
  printf("Results for %s \n--------------\n", name);
  printf("cf: %g+/-%g and chi2/ndf: %g\n", cfmin, cferror, chisqmin);
  
  char Opbase[150];
  sprintf(Opbase, "/home/sza5154/Research/feedersim/Theory/Output1024"); // this
  generate_theory_curves(dc, q, ftNL, Opbase, name, cf);
 
 return chisqmin;
 
}

void generate_theory_m3(double dc, char* Opbase) {
  FILE *fp;
  char fname[150];
  sprintf(fname, "%s/theorym3f.dat", Opbase);
  fp = fopen(fname, "w");
  double q, ftNL, nuc1, M, m3, rat;
  q=0.00003;
  int z;
  for (z=0; z<3; z++) {
    for (ftNL=10000000; ftNL<300000000; ftNL=ftNL*1.5) {
       for (M=1e13; M<1e16; M=M*1.1) {
       nuc1=nuc(dc, M, calA(q*ftNL), q, ftNL, z);
       m3 =  M3R(M, calA(q*ftNL), q, ftNL);
       rat =(1.0+0.66*m3)*ratio(M, dc, q, ftNL, z, 5);
       if (nuc1<5.2 && rat>1.0) {
        fprintf(fp, "%g %g %g\n", nuc1,m3, rat);
      }
  }
  }
}
  
  fclose(fp);
  sprintf(fname, "%s/theorym3h.dat", Opbase);
  fp = fopen(fname, "w");
  q = 1;
  for (ftNL=1.0; ftNL<1000; ftNL=ftNL*1.1) {
    nuc1=nuc(dc, M, calA(q*ftNL), q, ftNL, z);
    fprintf(fp, "%g %g %g\n", nuc1, M3R(M, calA(q*ftNL), q, ftNL), (1.0+0.29*M3R(M, calA(q*ftNL), q, ftNL))*ratio(M, dc, q, ftNL, z, 5));
  }
  
  fclose(fp);
     
}

void generate_theory_curves(double dc, double q, double ftNL, char* Opbase, char *name, double cf) {
  int z;
  FILE *fp;
  FILE *fpall;
  double M, nuc1;
  char fname[150], fnameall[150];
  sprintf(fnameall, "%s/%s/mfn/theory.dat", Opbase, name);
  fpall = fopen(fnameall, "w");
  int N;
  for (z=0; z<3; z++) {
    for (N=3; N<6; N++) {
    sprintf(fname, "%s/%s/mfn/theory%dz%d.dat", Opbase, name,N,z);
    fp = fopen(fname, "w");
    for (M=1e13; M<4e15; M=M*1.05) {
      fprintf(fp, "%g %g\n", M, ratio(M, dc, q, ftNL, z, N));
    }
    fclose(fp);
    }
    
    sprintf(fname, "%s/%s/mfn/theoryz%d.dat", Opbase, name, z);
    fp = fopen(fname, "w");
    for (M=1e13; M<4e15; M=M*1.05) {
      fprintf(fp, "%g %g\n", M, ratio(M, dc, q, ftNL, z, 5));
    }
    fclose(fp);  
    
    sprintf(fname, "%s/%s/mfn/theorycfz%d.dat", Opbase, name, z);
    fp = fopen(fname, "w");
    for (M=1e13; M<4e15; M=M*1.05) { 
      nuc1=nuc(dc, M, calA(q*ftNL), q, ftNL, z);
      if (cf!=0) {
        fprintf(fp, "%g %g\n", M, cf*ratio(M, dc, q, ftNL, z, 5));
        fprintf(fpall, "%g %g %g %g\n", nuc1, cf*ratio(M, dc, q, ftNL, z,5), M3R(M, calA(q*ftNL), q, ftNL), M3R(M, calA(q*ftNL),q,ftNL)*He3(nuc1)/6.0);
      } else {
        fprintf(fp, "%g %g\n", M, (f2hier(M3Rh(M,calA(q*ftNL),q, ftNL))*ratioh(M, dc, q, ftNL, z)+f2feed(M3Rf(M,calA(q*ftNL),q,ftNL))*ratiof(M, dc, q, ftNL, z))+fM(M, dc, q, ftNL,z));
        fprintf(fpall, "%g %g %g %g\n", nuc1, f2hier(M3Rh(M,calA(q*ftNL),q,ftNL))*ratioh(M, dc, q, ftNL, z)+f2feed(M3Rf(M, calA(q*ftNL),q, ftNL))*ratiof(M, dc, q, ftNL, z), M3R(M, calA(q*ftNL), q, ftNL), M3R(M, calA(q*ftNL), q, ftNL)*He3(nuc1)/6.0);
      }
    }
    fclose(fp);
    
    // also generate f1M curves
    
    sprintf(fname, "%s/%s/mfn/f1Mz%d.dat", Opbase, name, z);
    fp = fopen(fname, "w");
    for (M=1e13; M<4e15; M=M*1.05) {
      fprintf(fp, "%g %g\n", M, fM(M, dc, q, ftNL, z));
    }
    fclose(fp);
    
    mfn(dc, q, ftNL, name, 12., cf);
  }
  fclose(fpall);

}

double calA(double qftNL) {
  return 8.20092e-10*1.0/(0.5 + 6.325e-9*sqrt(6.25e15 + 3.2e8*pow(qftNL,2.0)));  // this factor has been checked to agree with the ratio measured directly from IC code; also, this factor 8.20092e-10 depends on the exact k0 used, but is irrevant for the mass function case, as what matters is the ratio that changes from the case of the factor that matches the sigma8, and this formula agrees with the ratio measured directly!
}

double mfn(double deltac, double q, double ftNL, char *name, double mu, double cf) {
//  printf("%s\n", name);
  char fname[150], string[500], fnameerror[150];
  FILE *fp;
  FILE *fperror;
  double M, rat, ratfNL, ratM3, A;
  A= calA(q*ftNL);
  int which;
  int ndf,z, ndfall=0;
  // open the relevant file for reading
  char Opbase[200];
  //sprintf(Opbase, "/home1/02539/sza5154/feedersim/Theory/Output1024");
  sprintf(Opbase, "/home/sza5154/Research/feedersim/Theory/Output1024"); 
  double chi2allz=0;
  double nuc1;
  FILE *fpall;
  char fnameall[150];
  sprintf(fnameall, "%s/%s/mfn/mbin.nuc.ratio.dat", Opbase, name);
  fpall=fopen(fnameall, "w");
  
  FILE *fnu;
  char fnuname[150];
  
  for (z=0; z<3; z++) {
    sprintf(fname, "%s/%s/mfn/mbinz%d.dat.ratio.avg", Opbase, name, z);
    fp = fopen(fname, "r");
    float col1, col2, col4, chi2=0.;
    ndf = 0;
    sprintf(fnameerror, "%s/%s/mfn/errorz%d.dat", Opbase, name, z);
    fperror = fopen(fnameerror, "w");
    
    sprintf(fnuname, "%s/%s/mfn/mbin.nucz%d.dat", Opbase, name, z);
    fnu = fopen(fnuname, "w");
    while (!feof(fp)) {
      if (fgets(string, 499, fp)) {
	sscanf(string, "%g %g %*g %g", &col1, &col2, &col4);
	M = col1*9.648e11; // in h^-1 M_sun
 	nuc1=nuc(deltac,M, A, q, ftNL,z); //printf("dc: %g, M: %g, A: %g, q: %g, ftNL: %g, nuc1: %g nuc1d: %g\n", deltac, M, A, q, ftNL,nuc1, nucD(deltac, M, A, q, ftNL, z));
 	fprintf(fpall,"%g %g %g %g %g\n", nuc1, M3R(M, calA(q*ftNL), q, ftNL), col2, col4, cf*ratio(M,deltac,q,ftNL,z*1.0,5));
//        Bf = 0.5+(nuc1-1.5)*3.0/4.0;
	if (col1>0 && nuc1 < mu) {
	  //if (M>2e14) { //try only the highest ones!
	  rat = cf*ratio(M,deltac,q,ftNL,z*1.0,5);
//	  printf("%g\n", rat);
	  chi2+=pow((rat-col2)/col4,2.0);
//	  if (ftNL==99) printf("dc: %g, z:%d, nuc: %g, ratth: %g, ratio3: %g, col2: %g\n", deltac,z, nuc1, rat, ratio3h(M,A,deltac, q, ftNL, z), col2);
 //if (ftNL==500) printf("dc: %g, z:%d, nuc: %g, ratth: %g, col2: %g\n", deltac,z, nuc1, rat, col2);
	  fprintf(fperror, "%g %g\n", nuc1, (col2-rat)/col2);
	  ndf++; //}
	  
	  fprintf(fnu, "%g %g %g %g %g %g\n", M, nuc1, M3R(M, calA(q*ftNL), q, ftNL), col2, col4, cf*ratio(M, deltac, q, ftNL, z*1.0, 5));
	  //printf("r3: %g, r:4 %g, r5: %g\n", ratio3(M, A, deltac, q, ftNL, z), ratio4(M, A, deltac, q, ftNL,z), ratio5(M, A, deltac, q,ftNL, z));
	}
      }
    }
   // printf("nz1: %d, \n", nz1);
    //printf("%s, z:%d, chi2: %g, ndf: %d, chi2/ndf: %g\n", name, z, chi2, ndf, chi2/ndf);
    fclose(fp);
    fclose(fperror);
    fclose(fnu);
    chi2allz+=chi2;
    ndfall+=ndf;
  }
  fclose(fpall);
  NDFall=ndfall;
  //printf("ndf: %d, nz1: %d \t ", ndfall, nz1);
  return chi2allz/(ndfall-1.0);
}

double f2hier(double M3) {
  return 1.0+0.29*M3;
}

double f2feed(double M3) {
  return 1.0+0.66*M3;
}

double fM(double M, double deltac, double q, double ftNL, double z) {
  double A;
  A = calA(q*ftNL);
  return (nucD(deltac, M, A, q, ftNL, z)/nucD(deltac, M, calA(0), 0, 0,z)) * exp(-0.5*(pow(nuc(deltac, M, A, q, ftNL, z), 2.0)-pow(nuc(deltac, M, calA(0), 0, 0, z), 2.0)));
}

double ratioh(double M, double deltac, double q, double ftNL, double z) {
  double A;
  A = calA(q*ftNL);
  return fM(M, deltac,q,ftNL,z)*(ratio3h(M,A,deltac,q,ftNL,z)+ratio4h(M,A,deltac,q,ftNL,z)+ratio5h(M,A,deltac,q,ftNL,z));
}

double ratiof(double M, double deltac, double q, double ftNL, double z) {
  double A;
  A = calA(q*ftNL);
  return fM(M, deltac,q,ftNL,z)*(ratio3f(M,A,deltac,q,ftNL,z)+ratio4f(M,A,deltac,q,ftNL,z)+ratio5f(M,A,deltac,q,ftNL,z));
}

double ratio(double M, double deltac, double q, double ftNL, double z, int N) {
  double A;
  A = calA(q*ftNL);
  if (N==3) {
    return fM(M, deltac, q, ftNL, z)*(1.0+ratio3h(M, A, deltac, q, ftNL, z)+ratio3f(M, A, deltac, q, ftNL, z));
  }
  else if (N==4) {
    return fM(M, deltac, q, ftNL, z)*(1.0+ratio3h(M, A, deltac, q, ftNL, z)+ratio3f(M, A, deltac, q, ftNL, z) + ratio4h(M,A,deltac,q,ftNL,z)+ratio4f(M,A,deltac,q,ftNL,z));
  }
  else {
    return fM(M, deltac, q, ftNL, z)*(1.0+ratio3h(M, A, deltac, q, ftNL, z)+ratio3f(M, A, deltac, q, ftNL, z) + ratio4h(M,A,deltac,q,ftNL,z)+ratio4f(M,A,deltac,q,ftNL,z) +ratio5h(M,A,deltac,q,ftNL,z)+ratio5f(M,A,deltac,q,ftNL,z));
}
}

double ratio3h(double M, double A, double deltac, double q, double ftNL, double z) {
  return M3Rh(M,A,q,ftNL)*He3(nuc(deltac, M, A, q, ftNL,z))/6.0-M3RDh(M,A,q,ftNL)*He2(nuc(deltac,M,A,q,ftNL,z))/(6*nucD(deltac,M,A,q,ftNL,z)) ;
}

double ratio3f(double M, double A, double deltac, double q, double ftNL, double z) {
  return M3Rf(M,A,q,ftNL)*He3(nuc(deltac, M, A, q, ftNL,z))/6.0-M3RDf(M,A,q,ftNL)*He2(nuc(deltac,M,A,q,ftNL,z))/(6.0*nucD(deltac,M,A,q,ftNL,z)); 
  
  //   return M3Rf(M,A,q,ftNL)*He3(nuc(deltac, M, A, q, ftNL,z))/6.0-M3RDf(M,A,q,ftNL)*He2(nuc(deltac,M,A,q,ftNL,z))/(6.0*nucD(deltac,M,A,q,ftNL,z))+ M2R(M, A, q, ftNL)*He2(nuc(deltac, M, calA(0), 0, 0, z)/2.0)-M2RD(M,A,q,ftNL)*He1(nuc(deltac,M,calA(0),0,0,z))/(2.0*nucD(deltac,M,calA(0),0,0,z)); 
}

double ratio4h(double M, double A, double deltac, double q, double ftNL, double z) {
  return  M4Rh(M,A,q,ftNL)*He4(nuc(deltac,M,A,q,ftNL,z))/24-M4RDh(M,A,q,ftNL)*He3(nuc(deltac,M,A,q,ftNL,z))/(24.0*nucD(deltac,M,A,q,ftNL,z))+pow(M3Rh(M,A,q,ftNL),2.0)*He6(nuc(deltac,M,A,q,ftNL,z))/72.0-M3Rh(M,A,q,ftNL)*M3RDh(M,A,q,ftNL)*He5(nuc(deltac,M,A,q,ftNL,z))/(36.0*nucD(deltac,M,A,q,ftNL,z));
}

double ratio4f(double M, double A, double deltac, double q, double ftNL, double z) {
  return  M4Rf(M,A,q,ftNL)*He4(nuc(deltac,M,A,q,ftNL,z))/24.0-M4RDf(M,A,q,ftNL)*He3(nuc(deltac,M,A,q,ftNL,z))/(24.0*nucD(deltac,M,A,q,ftNL,z));
}

double ratio5h(double M, double A, double deltac, double q, double ftNL, double z)
{
  return M5Rh(M, A, q, ftNL)*He5(nuc(deltac, M, A, q, ftNL, z))/120.0-M5RDh(M,A,q,ftNL)*He4(nuc(deltac,M,A,q,ftNL,z))/(120.0*nucD(deltac,M,A,q,ftNL,z))+M3Rh(M,A,q,ftNL)*M4Rh(M,A,q,ftNL)*He7(nuc(deltac,M,A,q,ftNL,z))/(24.0*6.0)+pow(M3Rh(M,A,q,ftNL)/6.0,3.0)*He9(nuc(deltac,M,A,q,ftNL,z))/6.0;
}

double ratio5f(double M, double A, double deltac, double q, double ftNL, double z)
{
  return M5Rf(M, A, q, ftNL)*He5(nuc(deltac, M, A, q, ftNL, z))/120.0-M5RDf(M,A,q,ftNL)*He4(nuc(deltac,M,A,q,ftNL,z))/(120.0*nucD(deltac,M,A,q,ftNL,z));
}

double nuc(double deltac, double M, double A, double q, double ftNL, double z) {
  double gf;
  if (z==0){gf = 0.76001;} else if (z==0.5) {gf = 0.59455;} else if (z==1) {gf = 0.473345;} else if (z==2) {gf = 0.327534;} else { gf=0.; printf("no growth function value!\n"); }
  return deltac*0.76001/(sqrt(d2R(M,A,q,ftNL))*gf);
}

double nucD(double deltac, double M, double A, double q, double ftNL, double z) {
  double gf;
  if (z==0){gf = 0.76001;} else if (z==0.5) {gf = 0.59455;} else if (z==1) {    gf = 0.473345;} else if (z==2) {gf = 0.327534;} else {  gf=0.; printf("no growth function value!\n"); }
  return (-0.380005*deltac*(A*I2R1D(M) +pow(A,2.)*pow(ftNL,2.)*pow(q,2.)*
  I2R2D(M)))/(gf*pow(A*I2R1(M) +pow(A,2.)*pow(ftNL,2.)*pow(q,2.)*I2R2(M),1.5));
}

/*double M2R(double M, double A, double q, double ftNL) {
  return d2R(M,A,q,ftNL)/d2R(M,8.20092e-10,0,0)-1;
}*/
double M3R(double M, double A, double q, double ftNL) {
  return M3Rh(M, A, q, ftNL)+M3Rf(M, A, q, ftNL);
}

double M3Rh(double M, double A, double q, double ftNL) {
  return d3R1(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),1.5);
}

double M3Rf(double M, double A, double q, double ftNL) {
  return d3R2(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),1.5);
}

double M3RDh(double M, double A, double q, double ftNL) {
  return (-1.5*pow(A,2)*ftNL*pow(q,2)*I3R1(M)*
      (A*I2R1D(M) + 
        pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*
       pow(q,2)*I2R2(M),2.5) + (pow(A,2)*ftNL*pow(q,2)*I3R1D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*
       pow(q,2)*I2R2(M),1.5);
}

double M3RDf(double M, double A, double q, double ftNL) {
  return (-1.5*pow(A,3)*pow(ftNL,3)*pow(q,3)*I3R2(M)*
      (A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),
     2.5) + (pow(A,3)*pow(ftNL,3)*pow(q,3)*I3R2D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),1.5);
}

// EXACT SCALING FROM MC CALCULATIONS
double M4Rh(double M, double A, double q, double ftNL) {
  return d4R1(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),2.0);
}

double M4Rf(double M, double A, double q, double ftNL) {
  return d4R2(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),2.0);
}

double M4RDh(double M, double A, double q, double ftNL) {
  return (-2*pow(A,3)*pow(ftNL,2)*pow(q,3)*I4R1(M)*
      (A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
    pow(A*I2R1(M) +  pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),3)
     + (pow(A,3)*pow(ftNL,2)*pow(q,3)*I4R1D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),2);
}

double M4RDf(double M, double A, double q, double ftNL) {
  return (-2*pow(A,4)*pow(ftNL,4)*pow(q,4)*I4R2(M)*
      (A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),3)
     + (pow(A,4)*pow(ftNL,4)*pow(q,4)*I4R2D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),2);
}

double M5Rh(double M, double A, double q, double ftNL) {
  return d5R1(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),2.5);
}

double M5Rf(double M, double A, double q, double ftNL) {
  return d5R2(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),2.5);
}

double M5RDh(double M, double A, double q, double ftNL) {
  return (-2.5*pow(A,4)*pow(ftNL,3)*pow(q,4)*I5R1(M)*
      (A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),
     3.5) + (pow(A,4)*pow(ftNL,3)*pow(q,4)*I5R1D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),
     2.5);
}

double M5RDf(double M, double A, double q, double ftNL) {
  return (-2.5*pow(A,5)*pow(ftNL,5)*pow(q,5)*I5R2(M)*
      (A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),
     3.5) + (pow(A,5)*pow(ftNL,5)*pow(q,5)*I5R2D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),
     2.5);
}


//NAIVE THEORY SCALING
/***
double M4Rh(double M, double A, double q, double ftNL) {
  return An(4)*pow(M3Rh(M,A,q,ftNL)/6.0,2.0);
}

double M4Rf(double M, double A, double q, double ftNL) {
  return Bn(4)*pow(M3Rf(M,A,q,ftNL)/8.0, 4.0/3.0);
}

double M4RDh(double M, double A, double q, double ftNL) {
  return An(4)*2*M3Rh(M,A,q,ftNL)*M3RDh(M,A,q,ftNL)/36.0;
}

double M4RDf(double M, double A, double q, double ftNL) {
  return Bn(4)*pow(M3Rf(M,A,q,ftNL), 1.0/3.0)*M3RDf(M,A,q,ftNL)/12.0;
}

double M5Rh(double M, double A, double q, double ftNL) {
  return An(5)*pow(M3Rh(M,A,q,ftNL)/6.0, 3.0);
}

double M5Rf(double M, double A, double q, double ftNL) {
  return Bn(5)* pow(M3Rf(M, A, q, ftNL)/8.0, 5.0/3.0);
}

double M5RDf(double M, double A, double q, double ftNL) {
 return Bn(5)*M3RDf(M,A,q,ftNL)*pow(M3Rf(M,A,q,ftNL)/8.0,2.0/3.0)/12.0;
}

double M5RDh(double M, double A, double q, double ftNL) {
  return An(5)*3.0*pow(M3Rh(M,A,q,ftNL)/6.0, 2.0)*M3RDh(M,A,q,ftNL)/6.0;
}



***/

/*double M2RD(double M, double A, double q, double ftNL) { // check this again!
  return -(((A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*
           I2R2(M))*I2R1D(M))/pow(I2R1(M),2)) + (A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M))/I2R1(M);
}*/


double An(double n) {return tgamma(n-1.0)*pow(2.0, n-3.);}
double Bn(double n) {return tgamma(n)*pow(2.,n-1.);}

double d3R(double M, double A, double q, double ftNL) {
  return d3R1(M, A, q, ftNL) + d3R2(M, A, q, ftNL);
}

double d3R1(double M,double A,double q,double ftNL){return q*q*ftNL*A*A*I3R1(M);}
double d3R2(double M,double A,double q,double ftNL){return pow(q*ftNL*A,3.0)*I3R2(M);}

double d2R(double M, double A, double q, double ftNL) {
  return d2R1(M, A, q, ftNL) + d2R2(M, A, q, ftNL);
}

double d2R1(double M,double A,double q,double ftNL){return A*I2R1(M);}
double d2R2(double M,double A,double q,double ftNL){return pow(q*ftNL*A,2.0)*I2R2(M);}
double d4R1(double M, double A, double q, double ftNL) {
  return pow(q*A, 3.0)*pow(ftNL, 2.0)*I4R1(M);
}
double d4R2(double  M, double A, double q, double ftNL) {
  return pow(q*ftNL*A, 4.0)*I4R2(M);
}
double d4R(double M, double A, double q, double ftNL) {
  return d4R1(M, A, q, ftNL)+d4R2(M, A, q, ftNL);
}
double d5R1(double M, double A, double q, double ftNL) {
  return pow(q*A, 4.0)*pow(ftNL, 3.0)*I5R1(M);
}
double d5R2(double M, double A, double q, double ftNL) {
  return pow(q*ftNL*A, 5.0)*I5R2(M);
}
double d5R(double M, double A, double q, double ftNL) {
  return d5R1(M, A, q, ftNL)+d5R2(M, A, q, ftNL);
}

double I2R1(double M) {
  return (exp(211.95077430721494 + 53867.440698140446/pow(log(M),2) - 5707.298843590207/log(M))/
   pow(M,2.6982879649363776))/8.20092e-10;
}

double I2R2(double M) {
  return (exp(217.45099975669766 + 59994.611272720154/pow(log(M),2) - 6345.953357871641/log(M))/
   pow(M,2.999736998687599))/pow(8.20092e-10,2.);
}

double I3R1(double M) {
  return (exp(317.19706381126184 + 82519.39939502515/pow(log(M),2) - 8742.644240185558/log(M))/pow(M,4.153235049094334))/pow(8.20092e-10,2.0);
}

double I3R2(double M) {
  return (exp(328.5711956827565 + 90330.74012138607/pow(log(M),2) - 9555.216663852534/log(M))/ pow(M,4.52560678339216))/pow(8.20092e-10,3.0);
}

double I4R1(double M) {
  return (exp(435.05201222490723 + 114656.39541961903/pow(log(M),2) - 12128.556880039461/log(M))/pow(M,5.728602564734952))/pow(8.20092e-10,3.0);
}

double I4R2(double M) {
  return (exp(439.2909589964454 + 120387.65009836892/pow(log(M),2) - 12735.992294456157/log(M))/pow(M,6.02833747413132))/pow(8.20092e-10, 4.0);
}

double I5R1(double M) {
  return (exp(538.3664770385856 + 142274.36306453458/pow(log(M),2) - 15070.183565030553/log(M))/pow(M,7.153622574218454))/pow(8.20092e-10,4.0);
}

double I5R2(double M) {
  return (exp(536.0900011367548 + 145848.2928865423/pow(log(M),2) - 15474.924210049874/log(M))/ pow(M,7.387097398898808))/pow(8.20092e-10,5.0);
}


double I2R1D(double M) {
  return 1.2193753871516855e9*((-2.6982879649363776*
        exp(211.95077430721494 + 53867.440698140446/pow(log(M),2) - 
          5707.298843590207/log(M)))/pow(M,3.6982879649363776) + 
     (exp(211.95077430721494 + 53867.440698140446/pow(log(M),2) - 
          5707.298843590207/log(M))*(-107734.88139628089/(M*pow(log(M),3)) + 
          5707.298843590207/(M*pow(log(M),2))))/pow(M,2.6982879649363));
}

double I2R2D(double M) {
  return 1.4868763347913226e18*((-2.999736998687599*
        exp(217.45099975669766 + 59994.611272720154/pow(log(M),2) - 
          6345.953357871641/log(M)))/pow(M,3.999736998687599) + 
     (exp(217.45099975669766 + 59994.611272720154/pow(log(M),2) - 
          6345.953357871641/log(M))*(-119989.22254544031/(M*pow(log(M),3)) + 
          6345.953357871641/(M*pow(log(M),2))))/pow(M,2.999736998687599));
}

double I3R1D(double M) {
  return 1.4868763347913226e18*((-4.153235049094334*
        exp(317.19706381126184 + 82519.39939502515/pow(log(M),2) - 
          8742.644240185558/log(M)))/pow(M,5.153235049094334) + 
     (exp(317.19706381126184 + 82519.39939502515/pow(log(M),2) - 8742.644240185558/log(M))*
        (-165038.7987900503/(M*pow(log(M),3)) + 8742.644240185558/(M*pow(log(M),2))))/
      pow(M,4.153235049094334));
}

double I3R2D(double M) {
  return 1.813060406382848e27*((-4.52560678339216*
        exp(328.5711956827565 + 90330.74012138607/pow(log(M),2) - 9555.216663852534/log(M))
        )/pow(M,5.52560678339216) + (exp(
         328.5711956827565 + 90330.74012138607/pow(log(M),2) - 9555.216663852534/log(M))*
        (-180661.48024277214/(M*pow(log(M),3)) + 9555.216663852534/(M*pow(log(M),2))))/
      pow(M,4.52560678339216));
}

double I4R1D(double M) {
  return 1.813060406382848e27*((-5.728602564734952*
        exp(435.05201222490723 + 114656.39541961903/pow(log(M),2) - 
          12128.556880039461/log(M)))/pow(M,6.728602564734952) + 
     (exp(435.05201222490723 + 114656.39541961903/pow(log(M),2.) - 
          12128.556880039461/log(M))*(-229312.79083923806/(M*pow(log(M),3)) + 
          12128.556880039461/(M*pow(log(M),2))))/pow(M,5.728602564734952));
}

double I4R2D(double M) {
  return 2.2108012349624772e36*((-6.02833747413132*
        exp(439.2909589964454 + 120387.65009836892/pow(log(M),2) - 
          12735.992294456157/log(M)))/pow(M,7.02833747413132) + 
     (exp(439.2909589964454 + 120387.65009836892/pow(log(M),2) - 
          12735.992294456157/log(M))*(-240775.30019673784/(M*pow(log(M),3)) + 
          12735.992294456157/(M*pow(log(M),2))))/pow(M,6.02833747413132));
}

double I5R1D(double M) {
  return 2.2108012349624772e36*((-7.153622574218454*
        exp(538.3664770385856 + 142274.36306453458/pow(log(M),2) - 
          15070.183565030553/log(M)))/pow(M,8.153622574218453) + 
     (exp(538.3664770385856 + 142274.36306453458/pow(log(M),2) - 
          15070.183565030553/log(M))*(-284548.72612906917/(M*pow(log(M),3)) + 
          15070.183565030553/(M*pow(log(M),2))))/pow(M,7.153622574218454));
}

double I5R2D(double M) {
  return 2.6957966117977952e45*((-7.387097398898808*
        exp(536.0900011367548 + 145848.2928865423/pow(log(M),2) - 
          15474.924210049874/log(M)))/pow(M,8.387097398898808) + 
     (exp(536.0900011367548 + 145848.2928865423/pow(log(M),2) - 15474.924210049874/log(M))*
        (-291696.5857730846/(M*pow(log(M),3)) + 15474.924210049874/(M*pow(log(M),2))))/
      pow(M,7.387097398898808));
}


// the Hermite polynomials

double He1(double x) {return x;}
double He2(double x) {return x*x-1.0;}
double He3(double x) {return pow(x,3.0)-3.0*x;}
double He4(double x) {return pow(x,4.0)-6.0*x*x+3.0;}
double He5(double x) {return pow(x,5.0)-10.*pow(x,3.)+15.*x;}
double He6(double x) {return pow(x,6.0)-15.*pow(x,4.)+45.*x*x-15.;}
double He7(double x) {return pow(x,7.0)-21.*pow(x,5.)+105.*pow(x,3.)-105.*x;}
double He8(double x) {return pow(x,8.0)-28.*pow(x,6.)+210.*pow(x,4.)-420.*pow(x,2.)+105.;}
double He9(double x) {return pow(x,9.0)-36*pow(x,7.0)+378.0*pow(x,5.0)-1260.0*pow(x,3.0)+945.0*x;}
